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ABSTRACT 


A method for determining a region of relative 
stability for a class of interpolation and extrapolation 
formulas is obtained. Some considerations are also 
given to the absolute stability of this class of formulas. 
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ON THE STABILITY OF NUMERICAL SOLUTIONS 


OF ORDINARY DIFFERENTIAL EQUATIONS 

By Robert N. Lea 
Manned Spacecraft Center 

SUMMARY 


A method for determining a region of relative stability for a class of interpola- 
tion and extrapolation formulas is obtained. Some considerations are also given to the 
absolute stability of this class of formulas. 


INTRODUCTION 


This paper contains a theoretical development of a procedure for determining a 
region of relative stability for a given difference equation. A definition of relative 
stability is given in a paper by Ralston (ref. 1) which is concerned with the problem of 
integrating a single ordinary differential equation. The following results are con- 
cerned primarily with the problem of relative stability for integrating systems of 
ordinary differential equations, although some considerations are given to absolute 
stability as defined by Dahlquist (ref. 2) and others. Dahlquist’s work, however, dealt 
with the existence of stable numerical solutions, whereas this procedure is concerned 
with restricting certain parameters in order to insure that the solution is stable. 


SYMBOLS 


A 

matrix of partials A = ^ a^ 

a i’ b i’ D p+l’ 6 ’ e 

real constants 

a ij 

8f i/3 y. 

C.,a,X 

complex constants 

c (n) 

class of functions with the property that their first n deriv- 
atives exist and are continuous 



E 


displacement operator 
characteristic polynomial 


F 

F r > F t ,p r , p t , o , cr t partial derivatives with respect to the subscript 


f,y 


f 

n 

h 


i, j, k, m, n, p 
P, Q> P, o 
R 


R XR 


S,v 


t 


x 

n 


vector functions 

components of f and y, respectively 

f {Vn) 

integration step- size 
integers 
polynomials 
real numbers 

ordered pairs of real numbers 
open sets in R x R 
real variable 
a + nh 


X. 

1 


an approximation to y^x n 

- >fn) 

characteristic root of A 


DEFINITIONS 

The following is concerned with the stability of the numerical solution of the 
system 

y* = f(x, y), y(a) = y Q (1) 
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obtained through the use of the difference equation 

Vrt + * • • • + = h ( b /n + k * b k-i r n*-l + ' ' ' + b 0 r n) 1 

where h is the integration step- size, x q = a + nh, y n is an approximation to y^x^' 
and F n = f(x n ,y n ). Associated with every equation of the form of equation (2) is a ; 
pair of polynomials 


p(E) = aQ + a^E + — 
o(E) = bQ + bjE + • * < 

Jr 

If E is considered to be the displacement operator E y^ = y n+ ^, then equation (2) 
can be represented by 


*k E 


+ b. E 
k 


(3) 


p(E)y n = hcr(E)f n 


(4) 


In the following, it is assumed that: 

(a) The coefficients a^, b. are real, and a^ / 0 (i = 1, . . . , 

(b) No common factor exists for p and a. 


for yeC 


Definition 1. The difference equation (4) is said to have order 
(p+2) 


k) 

p if, and only if, 


p(E)y(x) - ha(E)y'(x) = D p+1 y (p+1) (x)h p+1 + 0 (h p+2 ) 


Definition 2. The difference equation (4) is consistent if its order p satisfies 

p - 1 > 0. 
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Definition 3. The difference equation (4) is said to be stable (in the sense of 
Dahlquist (ref. 2)) if 

(a) The roots of p are located within or on the unit circle 

(b) The roots on the unit circle are distinct. 

Throughout the following, it is assumed that equation (4) satisfies assumptions (a) and 
(b) stated previously, and (a) and (b) of definition 3. The following theorems proved 
by Dahlquist (ref. 2) limit the order of stable methods of the type as equation (4). 

Theorem 1. A necessary and sufficient condition that y n converge to y(x) is 
that the difference equation (4) be stable and consistent. 

Theorem 2. The order of a stable operator cannot exceed k+2, where k is the 
degree of the polynomial p. Furthermore, if k is odd, the order cannot exceed k+1. 

That stability and consistency are necessary conditions for convergence can be 
verified by applying a difference equation that does not satisfy definition 3 to the initial 
value problem y’ = 0, y(0) = 0. (See' ref. 3. ) Stability and consistency are also suffi- 
cient in the sense that, if no errors other than truncation error were introduced, the 
approximate solution would converge to the true solution as h approaches zero. How- 
ever, since round-off and other errors are unavoidable, and since a nonzero step-size 
must be chosen, it is necessary to define the meaning of stability in a neighborhood of 
zero. 


Definition 4. Let A^, j = 1, . . . , m denote the characteristic roots of A = 
and a„ = 3f.^9y^, where f. and y^ are the ith and jth components of the vectors 
f and y, respectively. The difference equation (4) is stable if, for all A^ with nega- 
tive real parts, the roots of 

P(r) - ta(r) = 0, t = hA . (5) 

lie within or on the unit circle. 

The equation 

[p(E) - tff(E)]e n = 0 (6) 
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gives the error in the numerical solution where e n = Y n “ y^x^. Equation (5) is called 

the characteristic error equation. If the roots of equation (5) are distinct, the solution 
of equation (6) has the form 


n 


= c i r i 


n 


+c k r k 


n 


(7) 


where the r^ are roots of equation (5). Note that the solution to equation (6) would be 

different if equation (5) had multiple roots. However, the form of the solution will not 
be of any concern in the following discussion. 

In the solution of a differential equation obtained by solving an approximating 
difference equation, instability arises due to the presence of extraneous solutions, that 
is, solutions of the difference equation that are not related to the solution of the differ- 
ential equation. In reference 3, it is shown that, if h is chosen so that the roots of 
equation (5) satisfy the condition of definition 4, errors introduced at some stage 
of the integration procedure are damped out in subsequent stages. 

In certain types of problems, it is important that the errors introduced by extra- 
neous solutions damp out faster than the true errors in the approximating solution. In 
these instances, the relative stability of equation (4) is important. Before relative 
stability is defined, the principal solution must be defined. 

Definition 5. The continuous function r j(h), satisfying r ^(0) = 1 and equa- 
tion (5), will be called the principal solution. All other solutions are extraneous. 

The following definition is given by Ralston in reference 1 . 

Definition 6. A consistent numerical integration method of the type as equa- 
tion (4) is said to be relatively stable on the interval [a, b], which must include zero, if 
for all he[a,b], 


r^h) 


< 1; i = 2, . . . , k 


where r^ denotes the principal solution of equation (5), and r^, i = 2, . . ., k, 
denote the remaining k-1 solutions. Furthermore, when jr^ = |r.| then r . is 
simple. 
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LOCATION OF THE PRINCIPAL SOLUTION 


In this section the implicit function theorem for functions of two variables will 
be needed. For a given function F(t, r), F^. and F r will denote, respectively, 

3F/3t and 3F/3r. 

Theorem 3 (implicit function theorem). Let F(t, r) be a function defined on an 
open set S of R x R, and let (t^, r^JeS. 

Suppose that 

(1) F has continuous first partial derivatives on S, 

(2) F(t Q , r Q ) = 0, and 

<3) F r( t 0' r 0)^ 0 


Then there exists an open neighborhood V of t Q and a unique function reC^ 
satisfying the following conditions: 

If teV, then 

(1) F(t,r(t)) = 0 

(2) (t, r(t))eS 

(3) dr/dt = -F t/ /F r 

Consider the characteristic polynomial (eq.(5)), which shall be denoted by 

F(h, a, r) = 0 (8) 

where a is an arbitrary complex parameter. 

The requirement that equation (4) be consistent implies F(0, a, 1) = 0. Also, 
as a consequence of condition (b) of definition 3 imposed on the polynomial p, 

F r (0, a, 1) / 0. Hence, by theorem 3, there exists some open region V containing 

zero in which r^(h) is unique and differentiable. Hence, r ^(h) is the solution of the 

differential equation 

dr/dh = aa(r)j(p Y - hon^ (9) 

with initial conditions r(0) = 1, where and cr r denote 3p/3r and do/dr, respec- 
tively. 
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Theorem 4. The solution of equation (9) is uniquely defined through any point 
(h, r) except possibly at those points (h, r) where r is a solution of 


and 


P° r ~ °P r = 0 } 




( 10 ) 


Proof. Let Aiq, r^^ be a solution of equation (9). Then 


p(r 0 )- h o^( r o)=° (ID 

and, if the solution is not unique, then 

p r( r o)- V a r( r o ) =0 <12) 


Case 1. Assume and a ^ 

hence, Tq is a root of equation (10). 
p and cr, ¥ 0. 


( r o) 


P( r o) p r( r o) 

are not zero. Then -7—: = ) ( and, 

a ( r 0) a r( r 0) 

Now as a result of condition (b) placed on 


Case 2. Assume CT r ( r g) = 0- 

If ff r ( r Q) = 0) then P r ( r g) = ® and r g I s S HH a solution of equation (10). 

Since cannot be zero, then for any solution (h^, r^ of equation (8), 11 q 

must be defined by 


h 


°" ao ( r 6) 


(13) 


This completes the proof. 

Hence, the principal solution of the difference equation is given uniquely by 
equation (9) in any region that does not include a solution of equation (10). 
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DETERMINATION OF A REGION OF STABILITY 


The following well-known theorem is stated for completeness. 

Theorem 5 (Rouche). If f(z) and g(z) are analytic interior to a simple closed 
Jordan curve C, and if they are continuous on C, and |f(z)| < | g(z)| on C, then the 
function F(z) = f(z) + g(z) has the same number of zeros interior to C as does g(z). 

The following theorem is an immediate consequence of theorem 5. 

Theorem 6. Let C be a simple closed Jordan curve so that an nth degree poly- 
nomial P(z) has exactly n-1 roots in the interior of C and P(z) 4 0 on C. Fur- 
thermore, let Q(z) be a polynomial and X a real number so that 

|P(z)| > |XQ(z)| 


on C; then if a is any complex number such that |a| < |x|, then P(z) ± ojQ(z) 
has exactly n-1 roots interior to C. 

Now consider the polynomials p and o where the degree of p is n and the 
degree of o is less than or equal to n. 

Theorem 7. Let C be any simple closed Jordan curve having n-1 roots of 
p(r) in its interior, and the root r ^ = 1 on the exterior. Let X be chosen such that 

I P (r ) | > | Xc(r)| 


for all r on C. Then r^(h) (the principal solution) is always on or outside C for 
all a such that | a | | X | . 

Proof. Certainly there is a neighborhood of zero such that r^(h) is outside 

C since r^h) is a continuous function. Indeed, all the roots of p(r) - Xo(r) are 

continuous functions of the coefficients (see for example ref. 4). Now if r ^(h) is ever 

interior to C, then there must exist real numbers h' and e> 0 such that r ^(h') is 

on C, and either r^h* +6) ^or r^h' -6)) is interior to C for any 6 such that 

0< 6< e. Since the remaining n-1 roots are interior to C, there exists a 
neighborhood of each root contained completely within C. Now consider an 
arbitrary 6 such that r j(h* + 6) is in C. Then some other root ^(h* + 6) 

must be outside or on C which contradicts the fact that is continuous. 

The main result of this paper follows as a direct consequence of the previous 
theorem and is stated as corollary 7.1. 
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Corollary 7.1. If X is chosen as in theorem 7, and C of theorem 7 is taken 
as a circle with center at the origin and radius less than 1, then the approximate 
solution to the system (eq. (1)) obtained from equation (4) is relatively stable if h is 
chosen such that ha < |X| where a is a bound on the eigenvalues of the matrix of 


partial derivatives 



Another corollary concerning the absolute stability of the difference equation 
solution follows. 


Corollary 7. 2. Let C be any simple closed Jordan curve contained in the unit 
circle and excluding the point z=l and let X be chosen as in theorem 7. Then if 
I ho'l <|x| for all eigenvalues a of the matrix of partial derivatives, the absolute 

stability of the numerical solution is dependent only on the principal solutions of the 
system corresponding to each eigenvalue. 


AN EXAMPLE 


Before considering a simple illustration of the results of the preceding section, 
observe the following theorem. 

Theorem 8. If p(E) = hcr(E )f^ is a difference equation such that: 

(1) p has all zero roots with the exception of the principal root 

(2) r i? i = 1, . . . , n are the roots of a 

(3) 6 is a real number such that 0 < 6 < 1 

(4) a is a bound on the partial derivatives of A (the matrix of partial deriva- 

tives), 

then the difference equation is relatively stable if 



( 14 ) 


where h is the integration step-size. 
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In view of theorem 7, the proof of the above proposition is obvious since 
|p(z)| >|p(6)| and |o(z)| < n^|r.| + 6^ for Izl =6. 

Now consider the system of equations 


*i = f i ( x r x 2) 

X 2 = f 2( X l’ X 2) 




( 15 ) 


and let the difference equation to be used be the trapezoidal formula 

y n + l- y n = 5( f n + l +f n) < 16 > 

For the system (eq. (15)) the matrix A is 

0 1 ' 

_ -1 0 . 

whose eigenvalues are obviously bounded by a = 1. For the difference equation (16) 
the associated polynomials are 

p(z) = z - 1 


a(z) = ^ (z + 1) 


Consequently, by theorem 8, the numerical solution will be relatively stable if h< 1. 
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CONCLUSIONS 


The primary purpose of this paper has been to determine methods for finding 
stability regions for numerical solutions of systems of differential equations using a 
difference equation. Very good bounds on stability regions are known or can be com- 
puted with relative ease for single differential equations. Thus this special case has 
not been explicitly mentioned. Note, however, that by using the results of the paper an 
optimum region of relative stability can be obtained. A method of determining a region 
of relative stability for the general case is given . While in general this region will 
not be optimal, it is at least easy to compute. For any step-size h such that the 
principal solution of the difference equation is a good approximation to the true error 
equation, stability of the difference equation solution is insured if the proper conditions 
on the eigenvalues are satisfied. 


Manned Spacecraft Center 
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